% scripts to make input data set for laterally inhomogeneous
% lithosphere calcualtions for Neesha's GJI paper


%% Station locations
station = {'NWP' 41.10 159.96 5580 08/01/01  08/01/02 90 ;...
'T13' 24.98 139.30 4794 10/01/05  11/30/06 90 ;...
'T14' 22.00 139.50 4945 10/01/05  11/30/06 90 ;...
'T15' 29.00 141.32 4026 10/01/05  11/30/06 90 ;...
'T16' 32.52 143.96 5408 11/01/07  11/30/08 90 ;...
'T18' 27.14 147.17 5594 11/01/06  11/30/07 90 };


station_lat = [ 41.10  24.98 22.00 29.00 32.52 27.14 ];
station_long = [ 159.96 139.30 139.50  141.32 143.96  147.17 ];
station_names = ['NWP';'T13';'T14';'T15';'T16';'T18'];

%% Create a lithospheric slab under phillipine sea

% load the existing lithosphere (ocean/cont) data

data = load('/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_180_360_modeling/sh_1D_Shimizu_lateral_lith_3K_.3K/lithosphere.txt');
station_lat = [ 41.10  24.98 22.00 29.00 32.52 27.14 ];
station_long = [ 159.96 139.30 139.50  141.32 143.96  147.17 ];
station_names = ['NWP';'T13';'T14';'T15';'T16';'T18'];
correct_data = reshape(data',[360,180])';
imagesc([0.5:359.5],[-89.5:89.5],flipud(correct_data));
hold on;
plot(station_long, station_lat,'c.');

%% assign continental values to a slab under phillippine sea

correct_data(90-30:90-19,132:146) = 0.1;
data1 = reshape(correct_data',[180,360])';
%% save the data 
 fid = fopen('/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_180_360_modeling/sh_1D_Shimizu_lateral_lateral_lith/lithosphere_with_phillipine_slab.txt','wt');

for i = 1:360,
for j = 1:180,
fprintf(fid,'%3.1f ',data1(i,j));
end;
fprintf(fid,'\n');
end;

%% Plot the Lithospheric Boundary Map

data1 = load('/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_180_360_modeling/sh_1D_Shimizu_lateral_lateral_lith/lithosphere.txt');
correct_data1 = reshape(data1',[360,180])';
imagesc([0.5:359.5],[-89.5:89.5],flipud(correct_data1))
set(gca,'YDir','normal')
% station_lat = [ 41.10  24.98 22.00 29.00 32.52 27.14 ];
% station_long = [ 159.96 139.30 139.50  141.32 143.96  147.17 ];
% station_names = ['NWP';'T13';'T14';'T15';'T16';'T18'];
% 
% hold on;
% plot(station_long, station_lat,'c.');
% text(station_long, station_lat,station_names)
% text(station_long, station_lat,station_names,'color','c')
set(gca,'FontSize',16);
set(gca,'xtick',[],'ytick',[]);
text(60,50,'3000 \Omega-m','FontSize',16,'color','w');
text(180,-10,'300 \Omega-m','FontSize',16,'color','w');
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Lithospheric_Boundary_Map.eps','eps2');
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Lithospheric_Boundary_Map.fig','fig');

%title('Lithospheric Boundary Map');
%%
% Now read the predicted data
File_1 = '/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_180_360_modeling/sh_1D_Shimizu_lateral_lateral_lith/tides_O1_180_360.res';
fid = fopen(File_1, 'r', 'n');
LBlock = fread(fid, 1, 'long');
NO = fread(fid, 1, 'long');
Period = fread(fid, 1, 'double');
np = fread(fid, 1, 'long');
nt = fread(fid, 1, 'long');
fseek(fid, LBlock*4, 'bof');
Hpr = fread(fid, [nt np], 'float')*4*pi*100.;
Hpi = fread(fid, [nt np], 'float')*4*pi*100.;
Htr = fread(fid, [nt np], 'float')*4*pi*100.;
Hti = fread(fid, [nt np], 'float')*4*pi*100.;
Hrr = fread(fid, [nt np], 'float')*4*pi*100.;
Hri = fread(fid, [nt np], 'float')*4*pi*100.;
% now project the prediction to the main field
% create unit vector in the direction of main field
F = sqrt(Bx_360(:).^2+By_360(:).^2+Bz_360(:).^2);
main = [Bx_360(:)./F By_360(:)./F Bz_360(:)./F];
% find the scalar anomaly by projecting the predicted vector
% on to the direction of main field
% Bx = -Ht, Bz = -Hr
B_scalar_real = dot([-Htr(:) Hpr(:) -Hrr(:)]',main');
B_scalar_imag = dot([-Hti(:) Hpi(:) -Hri(:)]',main');
% reshape the matrix back
Sclar_Anly_Real =  reshape(B_scalar_real,[180,360]);
Sclar_Anly_Imag =  reshape(B_scalar_imag,[180,360]);

% Now read the predicted data with lateral lithospheric inhomogenities
File_1 = '/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_180_360_modeling/sh_1D_Shimizu_lateral_lateral_lith/tides_O1_180_360_3K_300.res';
fid = fopen(File_1, 'r', 'n');
LBlock = fread(fid, 1, 'long');
NO = fread(fid, 1, 'long');
Period = fread(fid, 1, 'double');
np = fread(fid, 1, 'long');
nt = fread(fid, 1, 'long');
fseek(fid, LBlock*4, 'bof');
Hpr = fread(fid, [nt np], 'float')*4*pi*100.;
Hpi = fread(fid, [nt np], 'float')*4*pi*100.;
Htr = fread(fid, [nt np], 'float')*4*pi*100.;
Hti = fread(fid, [nt np], 'float')*4*pi*100.;
Hrr = fread(fid, [nt np], 'float')*4*pi*100.;
Hri = fread(fid, [nt np], 'float')*4*pi*100.;
fclose all
% create unit vector in the direction of main field
F = sqrt(Bx_360(:).^2+By_360(:).^2+Bz_360(:).^2);
main = [Bx_360(:)./F By_360(:)./F Bz_360(:)./F];
% find the scalar anomaly by projecting the predicted vector
% on to the direction of main field
% Bx = -Ht, Bz = -Hr
B_scalar_real = dot([-Htr(:) Hpr(:) -Hrr(:)]',main');
B_scalar_imag = dot([-Hti(:) Hpi(:) -Hri(:)]',main');
% reshape the matrix back
Sclar_Anly_Real_lateral =  reshape(B_scalar_real,[180,360]);
Sclar_Anly_Imag_lateral =  reshape(B_scalar_imag,[180,360]);



% Now read the predicted data with lateral lithospheric inhomogenities and
% a slab in the Phillippines
File_1 = '/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_180_360_modeling/sh_1D_Shimizu_lateral_lateral_lith/tides_O1_180_360_3K_300_slab.res';
fid = fopen(File_1, 'r', 'n');
LBlock = fread(fid, 1, 'long');
NO = fread(fid, 1, 'long');
Period = fread(fid, 1, 'double');
np = fread(fid, 1, 'long');
nt = fread(fid, 1, 'long');
fseek(fid, LBlock*4, 'bof');
Hpr = fread(fid, [nt np], 'float')*4*pi*100.;
Hpi = fread(fid, [nt np], 'float')*4*pi*100.;
Htr = fread(fid, [nt np], 'float')*4*pi*100.;
Hti = fread(fid, [nt np], 'float')*4*pi*100.;
Hrr = fread(fid, [nt np], 'float')*4*pi*100.;
Hri = fread(fid, [nt np], 'float')*4*pi*100.;
% now project the prediction to the main field
% create unit vector in the direction of main field
F = sqrt(Bx_360(:).^2+By_360(:).^2+Bz_360(:).^2);
main = [Bx_360(:)./F By_360(:)./F Bz_360(:)./F];
% find the scalar anomaly by projecting the predicted vector
% on to the direction of main field
% Bx = -Ht, Bz = -Hr
B_scalar_real = dot([-Htr(:) Hpr(:) -Hrr(:)]',main');
B_scalar_imag = dot([-Hti(:) Hpi(:) -Hri(:)]',main');

% reshape the matrix back
Sclar_Anly_Real_lateral_slab =  reshape(B_scalar_real,[180,360]);
Sclar_Anly_Imag_lateral_slab =  reshape(B_scalar_imag,[180,360]);

%% Now plot the shit

% define axes

theta = [0.5 : 1. : 180.]';
phi   = [0.5 : 1. : 360.]';

%load boundary map
load coast

% plot the map
figure(1);


imagesc(phi, theta, Sclar_Anly_Real - Sclar_Anly_Real_lateral, [-0.2 0.2] );
hold on; plot([long; 360+long], 90-[lat; lat], '-k');
set(gca, 'FontSize' , 16);
colorbar;
set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);
set(gca,'xtick',[],'ytick',[]);
title 'O1 difference between models with uniform and laterally inhomogeneous lithosphere - real part ';
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_180_360_O1_Scalar_Diff_Lateral_Real.eps','eps');
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_180_360_O1_Scalar_Diff_Lateral_Real.fig','fig');

figure(2);

imagesc(phi, theta, Sclar_Anly_Imag - Sclar_Anly_Imag_lateral, [-0.2 0.2] );
hold on; plot([long; 360+long], 90-[lat; lat], '-k');
set(gca, 'FontSize' , 16);
colorbar;
set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);
set(gca,'xtick',[],'ytick',[]);
title 'O1 difference between models with uniform and laterally inhomogeneous lithosphere - imag part ';
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_180_360_O1_Scalar_Diff_Lateral_Imag.eps','eps');
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_180_360_O1_Scalar_Diff_Lateral_Imag.fig','fig');

figure(3);



imagesc(phi, theta, Sclar_Anly_Real - Sclar_Anly_Real_lateral_slab, [-0.2 0.2] );
hold on; plot([long; 360+long], 90-[lat; lat], '-k');
set(gca, 'FontSize' , 16);
colorbar;
set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);
set(gca,'xtick',[],'ytick',[]);
title 'O1 difference between models with uniform and laterally inhomogeneous lithosphere with slab - real part ';
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_180_360_O1_Scalar_Diff_Lateral_slab_Real.eps','eps');
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_180_360_O1_Scalar_Diff_Lateral_slab_Real.fig','fig');

figure(4);

imagesc(phi, theta, Sclar_Anly_Imag - Sclar_Anly_Imag_lateral_slab, [-0.2 0.2] );
hold on; plot([long; 360+long], 90-[lat; lat], '-k');
set(gca, 'FontSize' , 16);
colorbar;
set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);
set(gca,'xtick',[],'ytick',[]);
title 'O1 difference between models with uniform and laterally inhomogeneous lithosphere with slab - imag part ';
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_180_360_O1_Scalar_Diff_Lateral_slab_Imag.eps','eps');
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_180_360_O1_Scalar_Diff_Lateral_slab_Imag.fig','fig');


